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Abstract 



£3 We study a message passing model, applicable also to traffic prob- 

^ lems. The model is implemented in a discrete lattice, where particles 

move towards their destination, with fluctuations around the minimal 
distance path. A repulsive interaction between particles is introduced 
O hi order to avoid the appearance of traffic jam. We have studied the 

^ parameter space finding regions of fluid traffic, and saturated ones, 

being separated by abrupt changes. The improvement of the system 
performance is also explored, by the introduction of a non-constant 
potential acting on the particles. Finally, we deal with the behavior 
of the system when temporary failures in the transmission occurs. 
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1 Introduction 



The formulation of adequate models for describing message passing through 
a Computer Network [1], the traffic flow in a city [2], adsorption of molecules 
in a crystal [3], or the mechanism of High Temperature Superconductors [4], 
is getting increasing attention in the last years. 

All these problems share that they can be studied from the point of 
view of the Random Walk Theory. Nevertheless, the high number of degrees 
of freedom makes infeasible to extract relevant information from analytical 
calculations. The alternative is the use of numerical simulations. 

One can think in formulating models, as close to the reality as possible, 
including as many degrees of freedom as one is able to handle. The result- 
ing approaches are too complicate to study, not only in what concerns the 
description of relevant phenomena, but also because they lack of predictive 
power. 

An alternative approach is to go to simpler formulations, less realistic, 
but retaining the main features of the physical system. Such models are 
described by few parameters, being in this way easier to handle, and a global 
study of the parameter space is feasible, without losing physical intuition. 

Our hope is to obtain relevant results on complex systems by studying rel- 
atively simple models as is done in Statistical Mechanics (SM) where systems 
as complex as a real ferromagnetic lattice share many properties with models 
as simple as the Ising Model. In SM the relation between both systems is 
understood by means of such techniques as the Renormalization Group. 

Unfortunately we are far from proving that kind of relation for traffic 
models but the study of simple model is anyway necessary. 

When we deal with problems involving information flow, it is necessary to 
face with problems related with congestion processes. These systems undergo 
a transition from a situation of fluid traffic to another one characterized by 
the traffic jam. This change, can be related to a phase transition in SM, and 
the tools developed for statistical systems used. The first of all is to search 
for the variables governing this transition. 

The transmission in the lattice, is conditioned by the maximum number 
of particles that a site is able to manage, the interaction between particles, 
and the geometric constraints such as the coordination number (number of 
neighbors). In what concerns the interaction, as said before, we hope that 
only the general features will be relevant. 

In our model we implement movement by assigning to each particle a 
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final point. In real traffic, the particles are often stopped by obstacles in 
their way. These obstacles are simulated by limiting the number of particles 
that a site can buffer simultaneously. 

At this point, we return to the discussion of the problem from the point 
of view of SM. Let us focus on the situation where the density is high. It is 
intuitively clear that a particle should avoid zones of traffic jam by allowing 
it to surround obstacles. The system will be globally clearer, and we expect 
the throughput to be improved, in spite of the temporary misrouting of the 
particles. This idea is not only beneficial in situations of traffic jam, but also 
in general: allowing fluctuations around the minimal distance path improves 
the characteristics of the traffic. 

In terms of SM, the above discussion means that by heating the system, 
(T ^ 0), we will obtain better performances [5]. 

In what concerns the interaction between particles, up to now we have in- 
troduced a contact interaction, the existence of obstacles, as sites completely 
occupied. The consequences of the contact interaction can be seen as an 
infinite potential acting on saturated sites. 

An appropriate interaction term, is obtained by assigning an electrical 
charge to every particle. The repulsion will move the particles away from 
regions of high charge. 

Repulsive forces are usually considered in the case of directed polymers [6]. 

We have introduced a term to control the size of the thermal fluctuations, 
and another one to simulate a repulsive force. The parameters controlling 
these terms are the temperature and the charge respectively. We will see 
how, by tuning both parameters, is possible to improve the performance of 
the transmission. 

2 The Model 

We consider a two-dimensional lattice with coordination number 4, and peri- 
odic boundary conditions. The particles live in the lattice sites, labeled by 
n = (n , ni) and can move from a site to one of its 4 nearer neighbors in every 
time step. The maximum number of particles that a site can buffer will be 
denoted by B, and will be kept fixed along all the simulation at B = 5. We 
denote the occupation number of the site n by a(n). In this notation, the 
particles are prevented of moving to sites with a in) = B. In addition we 
consider an infinite queue at each site for the messages waiting to be put on 
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the network. 

By analogy with SM systems, we will work with the inverse of the tem- 
perature, /3. We denote the charge of the particles by k and the probability 
of particle injection to the lattice by p. 

The dynamic of the system is the following: 

1. We choose at random a lattice site. 

2. A particle is added, with fixed probability p, to the queue of the site, 
waiting for being introduced in the lattice. 

3. If the queue of the site is not empty, and a(n) < B, a new particle is 
introduced in the lattice, and an endpoint is assigned to it randomly. 

4. All particles in the considered site, try to move towards one of its 4 
neighbors. For a given particle sitting on n we must assign a probability 
of jumping to each of its neighbors. This probability is given by: 

P(±fj) = A^exp(±/5sign(n / (U - n^) - Ka{n^)) , (1) 

where means the \x coordinate of the site n, and is the endpoint of 
the particle considered, and N is the normalization constant. To choose 
between possible destinations we use a Heat Bath [7] algorithm. 

5. The movement is granted if the chosen site has a(n) < B. Otherwise 
the particle remains at its original site until a new movement attempt 
in the next iteration. 

6. The particle disappears when moving to its endpoint. 

In this way, the lattice sites simulate nodes of transmission and reception 
of messages coming up from the outside. 

The factor multiplying /?, is a potential term. It implies a constant force 
acting on the particle and driving it to its endpoint. The k term produces 
repulsion between particles sitting in nearer neighbors sites. Obviously there 
is a wide range of potentials that could be considered, in order to produce 
more effective forces, and partial improvements will be expected. 

The less is p, the weaker is the effect of the interaction between the 
particles. For fixed (/?, k) the flow is feasible only for those values of p being 
less than a certain threshold (j3, k) dependent. Above this threshold, the 
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density becomes too high and the transmission process is prevented. We say 
the system is saturated. The saturation mechanism begins with the appear- 
ance of saturated domains, whose size grows collapsing the whole lattice. 

Our purpose is to quantify this threshold density, as well as to describe 
the flow properties along the parameter space (/?, n,p). 

Such a model can be considered as a first step towards an abstract model- 
ing of packet /message telecommunication networks. In such networks, pack- 
ets are routed by each node according to routing tables dynamically main- 
tained by the network so as to minimize some cost function along the traject- 
ory of the packets [8]. In the present model, this cost function is simply the 
traveled distance, and, as the network state does not change during simula- 
tions, routing tables need not be dynamically calculated. A natural extension 
of this model would be to allow variable link lengths between the nodes and 
to compute routing tables according to the path lengths. The role of the tem- 
perature parameter is to quantify how firmly the network will try to stick to 
its policy of minimizing the cost along a trajectory: as will be shown below, 
implementing the routing according to the cost minimization in a rigid way 
(zero temperature), although the best naive choice from an individual user 
point of view, can be detrimental to the global behavior of the network and 
thus to its collective utility (the total throughput offered by the network). 

Congestion is an unavoidable phenomenon in uncontrolled packet net- 
works, since packets introduced at one node have no guarantee of finding the 
necessary resources (available buffer space) in the transit nodes. If uncon- 
trolled, congestion results in packet losses which, in the case of data trans- 
mission, are detected by upper layers protocols (TCP for instance in the case 
of IP networks) and corrected by retransmission requests: packet losses may 
thus trigger more packets to be introduced in the network therefore amplify- 
ing the congestion. Two approaches are possible to avoid this vicious circle: 
reducing the source emission rate when a congestion is detected (this is the 
TCP approach, and it relies on the sources to behave differently to the con- 
gestion indications) or trying to control the congestion inside the network by 
preventing packets to access overloaded areas (this approach does not rely 
on any source behavior: in the extreme congestion case, the source is simply 
denied any access to the network). This last option, congestion control inside 
the network, is presently a very active area of research in the networking com- 
munity [9]. To implement such a control, the overload information should 
be transmitted by a node to its neighbors at least. This is precisely what 
is modeled by the repulsive charge term of the present model: access to a 
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loaded node is discouraged, and access to an overloaded node is simply for- 
bidden. It should be noted that, in the model as it stands now, a node is 
instantaneously aware of the load state of its neighbors. A more complex 
model should take into account the latency introduced by the node-to-node 
transit time, leading to a description of the congestion control by some kind 
of retarded potential. It will be particularly important to take this latency 
aspect into account for the modeling of high speed networks [10] , but this is 
left for a future work. 

3 Numerical Simulation 

We have done Monte Carlo (MC) simulations in order to study the parameter 
space. We present results obtained in a L x L lattice with L = 32 and 
periodic boundary conditions. The computations have been carried out on 
Workstations. 

The starting configuration is obtained by generating at each site a particle 
with probability p, therefore, the initial density is around p x L 2 . The ran- 
dom number generator is based on the one described in [11]. The time step 
is identified with a MC iteration, that is the update of L 2 lattice sites as 
described in section 2. 

The temporal evolution of the system exhibits a transient regime, char- 
acterized by the instability in the observables (see next section). After this, 
the system falls in an extremely long-living metastable state, where the flow 
properties do not change significantly with the temporal evolution. We say 
the system has reached the asymptotic regime. 

The time the system spends in the transient regime depends on the para- 
meter space point. In this regime, in the non-saturated region, the particle 
density is initially low and grows with the particle injection. 

This time also grows near the parameter space points where an abrupt 
change in the properties is observed (e.g. near the threshold density), reach- 
ing in this case up to 2 x 10 4 . 

At each value of the parameters, we have performed typically 8 x 10 4 
MC iterations. For j3 and k € (0, 0.4) the transient regime takes around 
5000 iterations, while above 0.5 for both parameters, this time falls to 400 — 
600 iterations. We have also made the simulations starting from different 
configurations, allowing the system to evolve until 2 x 10 5 iterations. We 
have computed the errors by calculating the dispersion between the results 
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obtained starting from different configurations. 



4 Observables 

A correct description of the system is obtained from the measure of relevant 
observables. From their temporal evolution we will be able to know when the 
asymptotic regime is reached, or even if this regime will be or not saturated. 
From the averaged value of these magnitudes we will obtain a quantitative 
description of the flow process. 

For a given configuration, we define the occupation M as: 

m = (i AO !>(*)■ (2) 

i=l 

The statistical average over the configurations (labeled by j) is: 

N 

(M)=KmJLl/N)Y l M j • (3) 
./ i 

being N the number of averaged configurations. 

During the time interval tj — U = At , rif particles will reach their end- 
point. We define the Band Width (B w ) as the number of particles arriving 
to the endpoint per time unit. 

Bw(At) = ^ . (4) 

The statistical average for B w is obtained from its mean value over a number 
of time intervals Nt - - 

N t 

{B w )= lim (l/N T )J2B w (j) . (5) 

J = l 

We can obtain a good measure of the system performance, from the mean 
time employed by the particles for reaching their destination, T M . We com- 
pute Tm in the time interval At, by adding the individual time spent by each 
particle (delay time), divided by nf. 

T M (At) = (6) 
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In the same way, we define the statistical average as over N T intervals as: 



N T 

{T M )= lim (1/N t )J2TmU) 



(7) 



The occupation frequency of a certain occupation number, a(n), is defined 
as the number of times that the occupation a(n) appears at any lattice site. 
The study of this quantity gives an indication of the required buffers 
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Figure 1: Temporal evolution of M bellow the critical injection p c . 



5 Phase Diagram 

We examine the parameter space ((3,n,p) searching for regions where there 
are sharp changes in the temporal evolution. At each k) value, there is a 
p value denoted by p c , such that for p < p c the asymptotic regime presents 
an stationary flow, and for p > p c the asymptotic regime is saturated and no 
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flow is possible. We plot the temporal evolution of M bellow p c (Figure 1) 
and above p c (Figure 2) for some values of the parameters. 
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Figure 2: Temporal evolution of M above p c 



This change is similar to a phase transition. The temporal evolution leads 
the system to one or another phase depending on the parameter space point. 
Once in the asymptotic regime, the non-saturated phase exhibits dynamic 
equilibrium: the number of injected particles, equals the number of arriving 
ones. This feature is reflected by (B w ) = p x L 2 . In this phase, (T M ) is 
constant, as well as (M) which is always less than B. 

The parameter space is divided in two regions by the surface defined by 
(f3, K,p c ). Figure 3 shows two sections, for fixed j3 and k respectively. 

Above the surface, after the transient regime, (M) = B, (B w ) = and 
(T M ) diverges. 

It is possible to give a simple interpretation of this congestion phenomenon; 
when the load increases so does the probability that two particles residing on 
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Figure 3: Sections through the phase space cube. On the left side is plotted 
p c versus j3, on the right one, p c versus k 



two neighboring saturated sites want to exchange their sites, but of course 
they cannot since one of the particles should move first and has nowhere to 
go as no resource is available at the desired site. This mutual blocking then 
tends to propagate to the other neighboring sites and above a given load, a 
complete congestion can develop at the scale of the network. Such a mutual 
blocking can also arise on a cycle of particles, although this is less probable 
for independent particles. 

Such mutual blocking is well-known in most no loss routing schemes such 
as wormhole routing and is called deadlock [12]. Deadlock avoidance is a very 
active research area, specially for massively parallel system interconnections 
[13]. 
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dependence 



The thermal fluctuations move the particle away from its minimal distance 
path. The higher is 0, the less important are these fluctuations. 

As shown in Figure 3 (left side), there is an interval of values where 
the reduction of fluctuations has a positive influence on the throughput: p c 
rises, as well as (B w ), while (T M ) decreases, as is shown is Figure 4. 

When going to great values, the situation does not persist. The absence 
of thermal fluctuations damage the throughput because the particles are not 
able of surrounding obstacles. As consequence p c and (B w ) decrease. (T M ) 
does not appreciably change from its minimal value, corresponding to the 
infinite directionality one. 
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Figure 2 shows how the thermal fluctuations influence the saturation 
mechanism. If they are important, M slowly grows until reaching B. If 
they are not significant, a sharp jump appears in the temporal evolution of 
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M, between its value in the transient regime and B. We conclude from this 
that the saturated domains grow faster in absence of fluctuations, as could 
be expected from the earlier discussions. 

In particular, it is intuitive that a mutual deadlock will last for a longer 
time if thermal fluctuations are disallowed (particles in mutual deadlock will 
repeatedly attempt to exchange their sites) hence creating a larger local 
congestion area which can eventually evolve towards a global congestion. 
This was clearly observed in [5] where only one buffer was available per site 
and this behavior also appears in this model, indicating that deadlock is more 
a consequence of a no loss routing scheme (movement is granted if and only 
if the available resource is available at the low end) than a consequence of 
insufficient resources. 
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Figure 5: (j3, K opt ) line. The errors are of the size of the k step measured, 
(Ak = 0.5) 
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k dependence 

The dependence of p c on k presents two different regions (see figure 3). 

In the first one, p c rises when increasing k. For these k values, the in- 
clusion of the repulsion term helps the system to avoid jammed regions. We 
denote by K opt the value in which p c reaches the maximum value. 

In the second region, k > K opt , p c decreases with k. From K opt on, the 
repulsion is too strong and it moves the particles far away from their minimal 
paths. As a result, the collapse appears for smaller injection densities. 

In Figure 5 we represent the value of n opt for some j5 values. The more 
restricted to its minimal path is the particle, the stronger is the repulsion 
needed in order to avoid obstacles in the lattice. 

Therefore, there is a maximum value for the particle injection supported 
at each j5 value, which is reached for n in the neighborhood of n opt . 
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6 General Behavior and Optimization 



Let us focus on the behavior of the other significant observables for the traffic 
flow, when j5 and n have been tuned in order to obtain the maximal p c . 

Behavior of B w 

As we have already pointed out, B w in the asymptotic regime is proportional 
to p by a factor L 2 , as it may be checked in figure 6. So at each (j3, k) value, 
the maximum for (Bw) is reached at p = p c . 

Above p c , (Bw) remains constant during the transient regime, and sharply 
falls to zero after it. Bw is a good measure in order to monitorize when the 
system has reached the asymptotic regime. 
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Behavior of T M 

The greater is p, the stronger is the effect of the interaction over the particles, 
they more and more turn aside their minimal path increasing (Tm) (Figure 
6). In Figure 7, we observe how (Tm) slightly increases with k until K opt . 
From K opt on, (Tm) increases with a high slope. 

We conclude that the better performances for the maximal injection sup- 
ported, and for the bandwidth, are obtained on the line (/?, n pt)- The repul- 
sion term damages Tm, the greater is k, the more time spend the particles in 
reaching their destination. Anyway, it is only for k > K opt that (Tm) increases 
in a dramatic fashion. 

Figure 8 shows the particle distribution in the simulation for j3 = 0.2 and 
j3 = 4.0. The plot of er(n)freq[er(n)] exhibits a maximum around (M), and 
also reveals a wider distribution of particles for small (5 values. 
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Figure 8: Distribution of the occupation number er(n) times the frequency 
of this state versus a(n) after 5 x 10 5 MC iterations. Transient regime con- 
tributions has been discarded. 
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Figure 9: j3 dependence of p c for the distance dependent force (dashed line). 
Dotted line represents the values obtained for the constant force 



7 Improving the throughput: force depend- 
ing on the distance 

Although our implementation of the physical system has been inspired by the 
search for the simplest model accomplishing the desired features, certainly, 
by taking into account details not included up to now, partial improvements 
are to be expected. As an example, we develop a possible improvement. In 
particular, we deal with the quenching of thermal fluctuations when they are 
not useful anymore, that is, in the neighborhood of the endpoints. 

We have seen how the j3 term acts as a potential approaching the particles 
to their destination. In the primary implementation, we have used a con- 
stant force, and therefore the particles support thermal fluctuations of the 
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same strength, no matter the remaining distance to the endpoint. We can 
expect a globally clearer lattice, if the particles distant a few sites from their 
destination are prevented of fluctuating. We will see how the inclusion of this 
feature do not spoil the good properties of the throughput, and also that the 
general behavior of the system is not altered. 

A possible implementation to take this fact into account is obtained by the 
introduction of a distance dependence in the probability distribution, in such 
a way that the contributions of the fluctuations decrease when decreasing the 
distance to the endpoint. 
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Figure 10: < T M > versus (5 for fixed k = 1 in p f« p c for the distance 
dependent force (dashed line) and constant force (dotted line). 

A straightforward way of doing that is to include a dependence on the 
relative approaching to the endpoint: The probability distribution reads now: 



P(±/i) = A^exp(±/?sign(n / jU - n^) - Ka(n^j) 



(8) 
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Where r n is defined by: 



"(«o,«l) 



y/(ji f o 



n ) 2 + (nfi — ni) 5 



(9) 



In Figure 9 the evolution of p c with (5 is shown. A global throughput 
improvement is reflected by higher p c values. This effect is more remark- 
able when the size of the thermal fluctuations is important (small j3 values), 
corroborating our first intuition on the effect of fluctuations in the steps 
preceding the endpoint. 

In Figure 10 we plot the (5 dependence of T M - We observe a global 
decrease in this time for all (5 values. 
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Figure 11: Delay distribution for the whole particles after 4 x 10 4 MC it- 
erations. Graphics on the left (constant force) must be compared with the 
right side ones (distance dependent force). The distributions on the top are 
calculated in (j3 = 0.2, k = l,p = 0.08) on the bottom they correspond to 
(/3 = 4,K = l,p = 0.008). 



The improvements concerning the delay time, are not limited to a smaller 
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(T M ). Figure 11 shows the delay distribution for all the particles, compared 
with the delay distribution obtained with the constant potential. We see at 
= 0.2 how the dispersion of the distribution strongly decreases when using 
the distance dependent force. So, the particles arrive in more similar times, 
increasing the uniformity and the reliability of the traffic flow. 



8 Fault tolerance 

As we have pointed out, the lattice sites simulate the control nodes of an 
information flow system. Needless to say that the sites are an idealization of 
nodes in real systems, because the possibility of communication failures is not 
allowed. Real nodes are subject to external factors, as technical constraints 
or outside influences that often damage the communication ability between 
some nodes. These nodes are then temporarily out of order, and they cannot 
communicate information nor receive it from the other ones. It is then inter- 
esting to have an idea of how robust the system is when temporary problems 
in the transmission exist. 

To implement communication problems, we prevent the information ex- 
change with n randomly chosen sites, during an interval At = 10 MC it- 
erations in the simulation. After this time interval, these n sites are again 
allowed to communicate, and other n are break at random 1 . 

In general the measured values for the observables will be shifted by an 
amount depending on n and on the parameter region. For each p value, we 
define the shift in the observable O as is value, relative to the one obtained 
in the ideal system: 

We have study the influence on the relevant parameters of communication 
failures, for n = 10, 100 and 200 in = 0.2 and = 4.0. 

In figure 12 is shown the comparative evolution of (M) for n = 10 and 
n = 100. 

The first consequence of a bad transmission is that the system must sup- 
port globally higher occupations. We measure smaller A r p c values the lower 
is 0, because, as has been already shown, the lattice supports higher occu- 
pation numbers when the fluctuations are important. 



1 We work in this section with the probability distribution (1) 
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A r (M) is almost zero for both (3 values, with n = 10 (around 1% of nodes 
out of order), while with n = 100, 200 A r (M) is always large. 
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Figure 12: (M) versus p, for some values of the faulting percentage. Dotted 
lines represent the (M) values obtained for the ideal lattice. On the top the 
parameter space point is (/3 = 4, k = 1), on the bottom the plotted point is 
(0 = 0.2, k = 1) 

In Table I we give the results obtained. As could be expected (T M ) 
increases for all n values, even though, it is almost zero for n = 10. 



9 Summary and Outlook 

We have studied a model useful for describing the relevant processes occurring 
in flow information systems, with immediate applications to network message 
passing and traffic problems in general. 
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Shift 


n = 10 


n = 100 


n = 200 


j3 = 4.0 k = 1.0 


A r p c 


0.09 


0.40 


0.58 


A r (T M )(~p c ) 


-0.05 


-0.38 


-0.89 


A r {M)(^p c ) 


0.00 


0.37 


0.55 


f3 = 0.2 k = 1.0 


A r p c 


0.04 


0.24 


0.41 


A r (T M )(^p c ) 


-0.00 


-0.08 


-0.37 


A r {M){^p c ) 


0.00 


0.43 


0.50 



Table 1: Shifted values of the some relevant observables for some values of 
the number of failures n. 

The introduction of the parameters (3 and k as controllers of the system 
behavior allows us to go a step further from purely descriptive models, be- 
cause we are able of giving prescriptions to improve the performance of the 
flow process. 

Deeper studies, skipped in this first work are also possible. Concretely, an 
accurate study of the scaling with L of the relevant magnitudes, as well as a 
detailed description of how the saturation time behaves, could lead us to the 
definition of quantities analogous to critical exponents. Also non-equilibrium 
states could be studied in order to monitorize the parameters controlling the 
saturation process. 

The introduction of the n parameter has implied that the particles are able 
of avoid jammed regions. The study has been limited to contact interactions, 
because the particles only see the occupation of its nearer neighbors. By 
informing the particles about the occupation of wider surrounding regions, 
improvements in the throughput are to be expected. 
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